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Understanding the PSCz Galaxy Power Spectrum with 
N-body Simulations 
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ABSTRACT 

By comparing the PSCz galaxy power spectrum with the results of nested pure dark- 
matter N-body simulations, we try to understand how infrared-selected galaxies pop- 
ulate dark-matter haloes, paying special attention to the method of halo identification 
in the simulations. We thus test the hypothesis that baryonic physics negligibly affects 
the distribution of galaxies down to the smallest scales yet observed. We are successful 
in reproducing the PSCz power spectrum on scales < 40 h Mpc , near our resolution 
limit, by imposing a central density cutoff on simulated haloes, which gives a rough 
minimum mass and circular velocity of haloes in which PSCz galaxies formed. 

Key words: large-scale structure of the universe - galaxies: haloes - cosmology: 
theory - galaxies: infrared - galaxies: formation - methods: A^-body simulations. 



1 INTRODUCTION 

The distribution of ice in the ocean can be difficult to study 
when only the tips of icebergs can be observed. We can cat- 
alogue the positions and redshifts of galaxies, and can ob- 
tain glimpses of the intergalactic environment by observing 
the Lyman-alpha forest, but the dark matter, the compo- 
nent which plays the largest role in our current paradigm 
of structure formation, remains obscure. In this paper, we 
try to connect underlying dark-matter ice to the iceberg tips 
(galaxies) we can see. 

Only in the last few years have we claimed to have a suc- 
cessful cosmological model, the 'concordance' ACDM model. 
Its loose ends do seem tieable, but some areas remain largely 
mysterious: for example, the formation of galaxies within 
dark-matter haloes, and the resulting relationship between 
them at the present epoch. Here, we investigate this rela- 
tionship through their distributions, most succinctly quan- 
tified by the power spectrum, or its Fourier dual, the cor- 
relation function. Currently, the most extensive measure- 
ment of a power spectrum of observed galaxies, ranging 
over 4.5 decades of wavenumber, is by Hamilton & Tegmark 
(2002, hereafter HT). It was made from the PSCz (Point 
Source Catalogue Redshift) catalogue (Saunders et al. 2000) 
of galaxies observed with the IRAS infrared satellite. There 
will soon be a flood of galaxy clustering data, for example 
from the Sloan Digital Sky Survey (SDSS, York et al. 2000) 
and the Two-Degree Field (Lewis et al. 2002) survey. Early 
data (e.g. Zehavi et al. 2002, Percival et al. 2001) suggest 
that clustering properties vary with galaxy morphology, lu- 



minosity, and colour. Here we restrict ourselves to PSCz 
infrared galaxies, but with excellent optical data differenti- 
ated by colour, an approach such as ours will soon be able to 
say more about the types of galaxies which inhabit different 
sorts of dark-matter haloes. 

As with other measurements of galaxy power spectra, 
for example from the APM galaxy survey (Baugh 1996), HT 
found a roughly power law form. This is strikingly different 
from the dark matter power spectrum in the current concor- 
dance ACDM cosmological model. Figurc^shows the PSCz 
power spectrum along with linear and non-linear power spec- 
tra for the concordance ACDM model, and also dark matter 
power spectra from our simulations. The non-linear dark 
matter spectrum traces the linear spectrum at large scales, 
but at smaller, non-linear scales, it rises above it because 
waves on the scale of collapsing structures grow faster than 
waves in the linear regime. At even smaller scales, virializa- 
tion slows growth, producing a downward inflection. 

The galaxy and dark-matter power spectra thus ap- 
pear to be biased with respect to each other (i.e. they are 
different); their different shapes indicate that bias is scale- 
dependent. Numerous attempts have been made to under- 
stand this. The halo model of large-scale structure (reviewed 
by Cooray & Sheth, 2002) assumes the existence of small, 
bound objects (haloes) which are clustered according to the 
linear power spectrum. Galaxy clustering statistics such as 
the correlation function can then be calculated as the sum 
of two terms describing pairs of galaxies from the same and 
from different haloes. This can make use of a Halo Occu- 
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Figure 1. Comparison of the dark matter power spectra from 
the simulations (thin solid curves) to the linear power spectrum 
(dotted curve) and, for the sake of illustration, an estimate of the 
non-linear power spectrum (dashed curve) evolved from the linear 
power spectrum using the method of Peacock & Dodds (1996). 
The PSCz power spectrum also appears, surrounded by a grey 
band. The boundaries of the 'non-linear' range of wavenumber 
were placed by hand, approximately at the two inflections in the 
dark-matter power spectrum. 



pation Distribution, HOD (Benson 2001; White, Hernquist 
& Springel 2001; Berlind & Weinberg 2002; Berlind et al. 
2002), which describes the number of galaxies which inhabit 
a halo of a given mass. In this context, a halo is defined 
as a region at least 200 (typical of virialization) times more 
dense than the background. The approximations of the halo 
model permit analytic models for the bias between galaxies 
and dark matter (e.g. Seljak 2000; Sheth & Jain 2002). Re- 
cently, a slight but statistically significant deviation from a 
power law in the projected correlation function of a prelim- 
inary SDSS sample was successfully modelled using a HOD 
in the halo model (Zehavi et al. 2003). 

Semi-analytic galaxy formation models (e.g. White & 
Frenk 1991; Kauffmann et al. 1999; Benson et al. 2000; 
Somerville et al. 2001, Mathis et al. 2002) and hydrodynamic 
simulations (e.g. Katz, Hernquist & Weinberg 1992; Cen & 
Ostriker 2000; Dave et al. 2000; Pearce et al. 2001; White, 
Hernquist, & Springel 2001; Yoshikawa et al. 2001) have also 
been successfully applied to the problem of bias. While these 
approaches attempt to give the relationship between galax- 
ies and dark matter directly, some of the galaxy formation 
prescriptions in semi-analytic models can be rather ad hoc, 
and it is not clear that hydrodynamic simulations correctly 
treat every piece of relevant baryonic physics. 

In this paper, we take a different tack: we directly fit the 
PSCz power spectrum to dark-matter-only N-body simula- 
tions. We are thus seeing how far we can take the assump- 
tion that on intergalactic scales, baryonic physics negligi- 
bly affects the clustering of haloes which contain observed 
galaxies. Similar studies, without fitting to specific observa- 
tions, and using different halo-finding algorithms, have been 
undertaken by Kravtsov & Klypin (1999) and Coh'n et al. 
(1999). 



We also do not employ a HOD framework. Such a statis- 
tical placement of galaxies in haloes is useful in constructing 
an analytic model or when faced with poor resolution, and 
can aid intuition. Our approach is more direct: if a halo is 
defined not as a region above a certain overdensity, but as 
a region gravitationally bound to a significant maximum in 
the density field, a definition which would admit subhaloes 
without their parent haloes, then the number of galaxies de- 
tected by a redshift survey inside a halo is either zero or 
one. 

We pay particular attention to the halo-finding algo- 
rithm we use to go from a dark matter distribution to a 
distribution of haloes, and to the effects of resolution on 
the algorithm. This is the most nontrivial step in compar- 
ing the results of simulations to the observed galaxy power 
spectrum, so it is important to be careful. 



2 METHOD 

2.1 The Simulations 

The PSCz power spectrum spans 4.5 decades of wavenum- 
ber; to replicate this dynamic range in an N-body simula- 
tion with sufficient mass resolution would be unfeasible. For 
this reason, and also to test for resolution effects, we ran 
four manageably sized 256 3 -particle simulations, of comov- 
ing box size 32, 64, 128, and 256 h" 1 Mpc, with an adaptive 
P 3 M code (Bertschinger 1991). A simulation of box size less 
than about 32 /i -1 Mpc would miss significant tidal forces 
from large-scale fluctuations, and also could not form large 
clusters that appear with low number densities in nature. 
The values of cosmological parameters we used in the simu- 
lation were from the concordance model of Wang, Tegmark 
& Zaldarriaga (2002): VL m = 0.34, fi A = 0.66, h = 0.64, and 
n — 0.93. We calculated the transfer function of the initial 
conditions with the code of Eisenstein & Hu (1999), which 
returned a value of erg = 0.63, the rms fluctuation of mass 
in spheres of radius 8 /i" 1 Mpc. 

There were two resolution issues to consider: mass res- 
olution and spatial resolution. The mass resolution, i.e. the 
mass of a particle in the simulation, depends on the number 
of particles per unit volume. (See Table|5|for particle masses 
for each simulation.) So, the mass resolution in the four sim- 
ulations necessarily changes with box size. However, we de- 
cided to use the same spatial resolution (softening length) for 
all simulations: 10 h" 1 kpc, roughly the lowest scale probed 
in the PSCz power spectrum. 

Unlike many cosmological simulations, our softening 
length was fixed in physical, not comoving (Eulerian, not La- 
grangian), coordinates. This means that at early epochs, the 
comoving softening length was larger, at maximum about 
1/6 the mean interparticle separation in the 32 h" 1 Mpc 
simulation, and less by factors of two in the others. The 
first haloes to collapse at z ~ 10 turn out to be quite im- 
portant in determining the fine structure of haloes at z = 0. 
The first haloes contain only a few particles at z f» 10, and a 
small comoving smoothing length can make their relaxation 
times tiny, resulting in an overproduction of small, relaxed 
structures (Moore 2001; Binney & Knebe 2002). By having 
a fairly large comoving softening length at early times, we 
hope to have mitigated this problem. 
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Table 1. Choice of initial conditions. Both the PSCz galaxy density and the 'Simulation' dark matter density are normalized to one. 





PSCz 


Simulation 


PSCz 


Simulation 


Corner boxes 


1.49255 


1.11328 


0.83325 


0.50146 




0.49410 


0.76367 


0.31982 


0.47900 




1.45612 


1.20386 


2.35201 


1.99805 




1.60033 


1.63550 


0.39196 


0.83447 


Central box 


1.36050 


1.10596 







We wanted the four simulated regions to be as sim- 
ilar as possible under the constraints of periodic bound- 
ary conditions, enabling us to compare structures in the 
four simulations to each other. Bertschinger (2001), build- 
ing on the work of Pen (1997), developed a simple formal- 
ism to generate Gaussian random fields with multiple lev- 
els of resolution, which we used to nest the initial condi- 
tions. Effectively, this means that the phases of fluctuations 
matched in the centres of each set of initial conditions. An 
animation depicting the nesting of the boxes can be found 
at http : //casa. Colorado . edu/~neyrinck/nesthalf .mpg 
The zone of agreement between two simulations of box size 
b and 6/2 is a central cube of side length roughly b/4; a 
greater zone of agreement is not possible since the smaller 
box has to obey periodic boundary conditions. 

We also selected our initial conditions so that the cen- 
tral region common to the four simulations would have a 
structure similar to the Local Group and its environs. In 
this way, we hoped to replicate some features of the way the 
PSCz observations were made, looking out from the Milky 
Way, in a slightly overdense region on the outskirts of a 
modest-sized supercluster. The highest-resolution informa- 
tion in the simulations is from the smallest box, just as the 
galaxies closest to each other in the PSCz catalogue came 
from regions close to the Milky Way. 

To obtain these special initial conditions, we generated 
ten sets of initial conditions for a 256 h~ x Mpc simulation on 
a 256 3 mesh, and ran them to the present epoch using a (fast 
but low-resolution) PM code. We counted particles in each 
cell of side length 8 Mpc, the traditional scale of non- 
linearity, over which we can trust the results of the quick PM 
code. This gave a 32 3 grid of dark matter density estimates p 
at the present epoch from each set of initial conditions. The 
same was done for the PSCz catalogue, binning galaxies on 
a 4 3 grid of 8 h' 1 Mpc cells with the Milky Way at the 
centre. This resulted in a galaxy number density (npscz) 
grid of total side length 32 ft -1 Mpc, far enough to enclose 
the local supercluster. We then compared all cubes 4 cells 
on a side from the simulations to the galaxy density grid. 
The comparison was made by minimizing the sum of (p — 
npscz) 2 over nine 16-/i _1 Mpc cubes: eight from dividing 
the 32 ft -1 Mpc cube into octants, and one more in the 
centre. We then shifted the best-fitting region to the centre 
of the 256 ft - 1 Mpc set of initial conditions before calculating 
the lower-box size initial conditions. Table Q shows these 
densities for the best fit. Unfortunately, this procedure did 
not result in structures with obvious visible similarity to the 
Local Group, but the statistical similarity is reassuring. 

Although similar, the inner regions of the four simula- 
tions were not identical after evolving them to the present 
epoch. Large-scale power caused bulk motion in the cen- 
tral region, moving it away and distorting it slightly from 
its original position. To assess this effect, we approximated 



it to first order with the sum of a translation and a lin- 
ear transformation, an asymmetric tensor which in general 
could include shear and rotation. We did not use this trans- 
formation in any analysis except in evaluating the similarity 
of the four simulations. The transformation can be written 
as 

rl 2) =C ijr V+ Si , (1) 

where is the ith coordinate in the central region of sim- 
ulation m, Cij is a deformation tensor, and Si is a translation 
vector. 

Figure shows particles initially in the central 16 h 
Mpc box from all simulations in the present epoch. The re- 
gion from the 32 h~ l Mpc simulation is untransformed, with 
the other three regions translated and deformed to fit best 
on to it. To calculate the shift and deformation, we compared 
the present-epoch positions of particles which had occupied 
the same place in the initial conditions; i.e., with the same 
Lagrangian positions. Since mass resolution changes across 
simulations, it was necessary to average together positions 
of particles (in sets of 2 3 , 4 3 , or 8 3 ) to compare positions in 
a smaller box size simulation to those in a larger one. 

We calculated the translation vector Si with {rf^ — r^), 
where the angular brackets denote an average over particles. 
As for the deformation tensor Cij, consider the quantity 
{( r i 2 ^ ~ s i) r k )■ Assuming that Eqn. Q holds, this equals 
Cijir^r^}. Thus, where B ik = <(rf > - a t )r^) and A jk = 

Cij = B ik (A x )kj- (2) 

The translation vectors Si from the 64, 128, and 256 
hT 1 Mpc simulations to the 32 had magnitudes 2.78, 4.43, 
and 5.92 h~ l Mpc, respectively. The deformation tensors 
dj were close to identity matrices, as one would hope; the 
diagonal elements were all between 0.89 and 1.06 except 
for one outlier at 0.81, and the off-diagonal elements had 
magnitudes below (mostly well below) 0.06. With the first- 
order adjustment, the agreement is still not perfect, but we 
expect the shear from large-scale waves outside the inner box 
to change slightly over its length; moreover, slight differences 
are likely to amplify when non-linear ly evolved. 

2.2 Halo finding 

To compare to the observed galaxy power spectrum, it is 
necessary first to find haloes in the set of dark-matter par- 
ticles returned by the simulations. Although they do rea- 
sonable jobs, no halo-finding algorithm (HFA) is perfect. In 
some analyses of N-body simulations, surprisingly little dis- 
cussion is given of the choice of HFA. 

The first step of most HFAs is to estimate the den- 
sity, a quantity which is not obviously defined given only 
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Figure 2. The region originally in the central 16 h~ 1 Mpc region of each simulation. Particles from the 32, 64, 128, and 256 h~ 1 
Mpc simulations appear in magenta, orange, green, and blue, respectively. The sets of particles from simulations of box size greater 
than 32 h" 1 Mpc have undergone a translation and linear transformation to lie on the set from the (untransformed) 32 h~ 1 Mpc 
simulation. All particles from the central region are shown for the 128 h~ 1 Mpc simulation, and the number densities of particles 
from other simulations are adjusted by factors of eight to match. Thus, 1/8 and 1/64 of the particles from the 64 and 32 Mpc 
simulations appear. In the 256 h~ 1 Mpc case, the number density was octupled by adding new particles which bisect lines joining 
particles adjacent on the initial mesh. The predominantly blue character of the large haloes at the bottom arises because the halo in 
the 256 h~ 1 Mpc simulation is simply in front of the others. The figure was produced using Nick Gnedin's IFRIT visualization tool, at 
http: //casa. Colorado . edu/~gnedin/ IFRIT/ 



a set of particles. DENMAX (Bertschinger & Gelb 1991) 
uses an Eulerian approach, calculating the density on a fine 
mesh by smoothing each particle with a Gaussian of a fixed 
size, called the smoothing length, r- smoo , on an effectively in- 
finitely fine mesh. A Lagrangian approach (HOP, Eisenstein 
& Hut 1998) uses a fixed number A^ens of nearest-neighbor 
particles to estimate the density at the position of each par- 
ticle, and also uses a few other parameters. DENMAX has 
a fixed spatial resolution, while HOP effectively has a fixed 
mass resolution. The results of both methods are strongly 
dependent on their free parameters, r smoo or N^ ells . 



Although DENMAX takes much more time to run than 
HOP, we ended up using a variant of DENMAX. We found 
that DENMAX is capable of finding smaller haloes than 
HOP, down to about ten particles. DENMAX works by mov- 
ing particles along density gradients until they are at a max- 
imum. It then uses a 'Friends-of-Friends' algorithm, finding 
clusters of moved particles closer than a small linking length 
(1/1024 times the box size) to each other. The last step is 
to 'unbind' iteratively any particles whose energies exceed 
the escape energy from their haloes. The output of DEN- 
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Figure 3. The effect of added DENMAX 2 haloes on the corre- 
lation function in the 32 h _1 Mpc simulation. The dashed curve 
is the PSCz correlation function, and the grey band shows its er- 
rors. The dotted curve is the best fit using just DENMAX haloes, 
and the solid curve is the best fit including DENMAX 2 haloes. 
The best-fitting central density cutoff was different in each case: 
for DENMAX, p c , min was 225 particles; for DENMAX 2 , it was 
122 particles. The arrows show the DENMAX and DENMAX 2 
smoothing lengths: 0.025 and 0.0125 /i -1 Mpc, respectively. 



MAX is a list of haloes with their masses (number of bound 
particles), and their position and velocity centroids. 

In DENMAX, a large smoothing length smears out close 
pairs of haloes, while a small smoothing length, with its 
higher density threshold, fails to include the outskirts of 
haloes in their mass, and also misses isolated, less-dense 
haloes. Gelb & Bertschinger (1994) discuss some effects of 
DENMAX resolution. Empirical tests have indicated that 
setting r smoo = 1/5, in units of the mean interparticle sep- 
aration, yields a halo mass spectrum similar to that given 
by the Press-Schechter (1974) formalism, a useful, though 
not omniscient, guide. This choice of r smoo makes some 
sense theoretically, too, since the spherical collapse model 
(Gunn & Gott 1972) predicts that a virialized object is 
5 w 180 times denser than the background in a standard fiat 
(p, m = 1) cosmology, and somewhat higher than that in a 
ACDM cosmology. Fiducially, regions of overdensity 200 and 
above are virialized, corresponding to a smoothing length of 
1/2003, about 1/5. 

We applied DENMAX with the canonical smoothing 
length to the results of each simulation, and calculated halo 
power spectra. The halo-finding resolution we obtained was 
rather poor; there was a small-scale downturn in each power 
spectrum starting at significantly larger scales than the sim- 
ulation's softening length. We therefore tried halving the 
DENMAX smoothing length to rsmoo = 1/10, which suc- 
ceeded in extending the power law in the correlation func- 
tion to smaller scales by about a factor of two. The smaller 
smoothing length evidently picked out subhaloes which the 
canonical smoothing length had merged together. A small 
smoothing length is desirable in detecting subhaloes within 
a halo, since in a halo, the spatial scales involved are smaller, 
and the background density is higher. 



We wanted to use higher resolution in higher-density 
regions without forsaking the advantages of the canonical 
r S moo in lower-density regions. We therefore used an algo- 
rithm which we call DENMAX 2 , in which DENMAX is run 
as normal with the canonical r smoo = 1/5, but then is ap- 
plied to each returned halo separately with r smoo = 1/10. 
Although this factor of two is rather arbitrary, it put the 
DENMAX 2 r smoo close to, but still above, the softening 
length in the 32 h~ 1 Mpc simulation, and it was a convenient 
choice with so many other factors of two being used, e.g. be- 
tween the box sizes. We could have kept the DENMAX 2 
r smoo the same in physical units between the boxsizes, but 
instead fixed it in interparticle units. This was for a cou- 
ple of reasons: we wanted to hold constant the ratio of the 
highest and lowest reliably measured wavenumber; also, in 
the larger-boxsize simulations, a fixed smoothing length in 
physical units would be tiny in interparticle units, and would 
thus encounter undesirably high Poisson noise. 

Figure 01 shows the effect of extra DENMAX 2 substruc- 
ture on the best-fitting (defined below) halo correlation func- 
tions in the 32 Mpc simulation; it extends the power law 
to scales about half as large, as one would expect from the 
halving of r amoo . The results are promising, but the choices 
of r smoo for DENMAX and DENMAX 2 are still arbitrary, 
indicating the desirability of a HFA without such free pa- 
rameters. 

From the list of DENMAX and DENMAX 2 haloes, we 
had to pick a subset which we thought could represent PSCz 
galaxies. Since DENMAX returns the mass of each halo, this 
was an obvious property to use to characterize haloes. How- 
ever, the mass of the largest DENMAX 2 subhalo is necessar- 
ily less than that of its parent DENMAX halo, both because 
it contains a subset of the parent halo's particles, and be- 
cause DENMAX 2 's smaller r rsmoo detects smaller structures. 
We thus could not compare DENMAX mass and DENMAX 2 
mass directly. We could have tried to 'correct' DENMAX 2 
masses to DENMAX levels, which are probably more physi- 
cally meaningful, since they approximate the virial mass. For 
example, we could have compared the DENMAX masses of 
each halo split by DENMAX 2 to the DENMAX 2 mass of its 
largest subhalo. However, it was occasionally unclear which 
subhalo represented the original halo. Also, it is not obvious 
that the relationship between DENMAX and DENMAX 2 
masses would be the same for satellite subhaloes and main 
subhaloes. 

We felt it was both simpler and more consistent to fit 
using another quantifier. Figure [I] shows that the central 
density of a halo, estimated by counting the number of par- 
ticles within a fixed radius r p of the halo's centre of mass as 
returned by DENMAX, is well-correlated to its DENMAX 
mass. Also, central density doesnot unfairly handicap dense 
but (as detected by DENMAX) unmassive haloes in their 
quest for inclusion in the halo list. For the 32 h" 1 Mpc sim- 
ulation, we used r p — 20 h~ x kpc. At twice the softening 
length, this was about the smallest scale at which we could 
expect to obtain a meaningful density estimate. Unfortu- 
nately, we could not use the same r p for the larger simula- 
tions, because the mass resolution became too poor, causing 
the density estimate within a small region to be dominated 
by Poisson noise. So, we simply increased r p in proportion 
to the box size of the simulation. 

Another quantity we could have used to characterize 
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Figure 4. Scatter plots of DENMAX mass (number of particles 
in the halo) versus central density p c (the number of particles 
within radius r p ) for all simulations. Haloes with p c = are 
shown with p c = 0.2. No haloes had a DENMAX mass under 
10 because DENMAX discards haloes with less than 10 particles. 
The solid grey line passing through the distribution shows the 
average mass at each p c . To obtain the mass at a particular p c ,0i 
we averaged together the masses of haloes with p c within 1/10 
of a dex of p c ,o- The thin black line shows the identity function, 
y = x, which would result if all, and only, particles counted in the 
DENMAX mass were within r D of the halo centre. 



haloes is maximum circular velocity v c . Figure shows a 
comparison of central density to v c for DENMAX haloes; 
the relationship is rather tight. To calculate v c , we used 
the procedure employed by the Bound Density Maxima 
HFA (Klypin et al. 1997). This procedure, described at 
http : //astro .nmsu. edu/~aklypin/PM/pmcode/node6 .html 
yields v c — -J GM(r) /r, evaluated at the radius of the first 
concentric shell about the halo where the overdensity inside 
dips below 200, and also checks its spherically averaged 
density profile to thwart structures hoping to inflate v c for 
nearby haloes. The maximum circular velocity is a useful 
quantity because it can seemingly be related to observed 
galactic circular velocities. However, it is fallible, relying as 
it does on spher averaged density profiles. In any case, the 
relationship between the two quantities is quite tight, so 
they are roughly interchangeable. 

We thus picked subsets of the halo list according to a 
lower density cutoff p c , m in', it is physically reasonable that 
the central density in a halo must exceed some threshold 
to house an observed galaxy. We first imposed the den- 
sity cutoff on the list of DENMAX haloes. If application 
of DENMAX 2 split a halo in the resulting subset into sub- 
haloes, we imposed the density cutoff on each of the sub- 
haloes. If none of these subhaloes exceeded the density cut- 
off, the original halo stayed in the list; otherwise, any dense- 
enough subhaloes (which almost always include the original) 
replaced the original halo in the list. 

We then calculated the haloes' power spectrum by bin- 
ning halo pairs by their separation, and then submitting the 
resulting correlation function to an FFT (actually FFTLog, 
Hamilton 2000). We also tried two other methods. In the 



Figure 5. Scatter plots of central density p c versus maximum 
circular velocity v c for all DENMAX haloes. The solid grey line 
passing through the distribution shows the average v c at each p c ■ 
To obtain v c at a particular p c ,0) we averaged together v c 's of 
haloes with p c within 1/10 of a dex of p c ,o- 



first, we calculated the density on a mesh, using the Nearest 
Grid Point interpolation scheme, and found the power spec- 
trum using a 3D FFT. To get to the smallest scales, Klypin 
(private communication, 2001) pointed out that if one di- 
vides a box of particles into octants (or some other number 
of equal parts) , and overlays all of the octants on each other, 
periodic boundary conditions will still be satisfied, and the 
power spectrum of the condensed box should be the same as 
that of the larger box. This approach did work rather well, 
but there were small discrepancies between power spectra 
from different octant overlayings, and it was not obvious 
how to combine them. Another method we tried used an 
unequally spaced FFT (Beylkin 1995), which uses multires- 
olution analysis (wavelets) to calculate the exact FFT of a 
set of delta functions in mass (i.e. particles), but a suffi- 
ciently large unequally spaced FFT required more memory 
than was convenient. These other methods agreed well with 
the correlation function technique, but we ended up using 
the correlation function technique because it is possible to 
calculate the exact correlation function of a relatively small 
number of haloes quickly with no resolution limit. Also, this 
technique is not subject to the vagaries of window functions 
which exist for standard FFTs. 

While it is more direct to calculate a correlation func- 
tion than a power spectrum from a simulation, the opposite 
is true for a redshift survey such as PSCz. This is because the 
power spectrum in directions transverse to the line of sight in 
a redshift survey is unaffected by redshift distortions, which 
is not true for the correlation function. So, in comparing sim- 
ulations to observations, it was always necessary to translate 
one set of data into the same space (either real or Fourier) as 
the other. Empirically, the correlation function varied more 
dramatically on the small-scale end with the density cutoff 
pc.min, and also it is easier to interpret directly in terms 
of physical pairs of haloes, so we used the correlation func- 
tion for fitting. It would have been better in principle to 
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Figure 6. The best-fitting power spectra (see Fig, 1131 for error 
bars) for all simulations. The PSCz power spectrum is the dashed 
curve with a grey error band; the best fits from the (from left to 
right) 256, 128, 64, and 32 h~ 1 Mpc simulations appear as solid 
curves. Triangles connected with dotted lines denote negative val- 
ues. 
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Figure 7. The best-fitting power spectra (error bars are similar 
to those in Fig. El for tnc 128, 64, and 32 h _1 Mpc simulations, 
constrained so that the halo number density matches. The PSCz 
power spectrum is the dashed curve with a grey error band. Tri- 
angles connected with dotted lines denote negative values. 



use the errors in the power spectrum, since they are more 
directly measured from PSCz. However, HT have not at 
present found a positive-definite covariance matrix for the 
PSCz power spectrum. Thus any comparison we make to 
simulations would not be completely rigorous anyway. We 
were still able to estimate the goodness of fit by ignoring 
cross-correlations among data points (just using HT's error 
bars), making a 'pseudo-x 2 ' (x 2 ) statistic. Where £ denotes 
the correlation function, and 6 is the box size of a simulation, 
we included f(r)'s with r between 6/256 and 6/2; with our 
bins ri varying as they did by a factor of \/2, this included 14 

points in the fit. We calculated x 2 = ( l(r ' o , ) ~ lpsCz(ri) ' 



3 "pso«(n) 

where £pscz(r) is the PSCz correlation function at r, and 

the average of 



2 



o"PSCz( r ) is the error in £pscz(r) = 
the upper and lower error bars as reported by HT. £pscz, 
and were logarithmically (or linearly, if adjacent data 
points straddled zero) interpolated if necessary. 



3 RESULTS 

Figure |S| shows the best individually- fit power spectra from 
each simulation. The resulting bias factors (see next para- 
graph) appear in Fig.|5] and the x 2 curves for the fits appear 
in Fig. Q3 Figure [TUl shows the x 2 curves as a function of halo 
number density. Fi gures ITI and [Til show an alternative collec- 
tive fit, still varying a central density cutoff, but constrain- 
ing the halo number density to be the same in the smallest 
three simulations; see H3.1.1l for further discussion. With the 
first of these fits, we have reproduced the form of the PSCz 
power spectrum on scales k < 40 h Mpc" 1 , or r > 0.05k" 1 
Mpc. The softening length was O.Ol/i -1 Mpc, so this is not 
much worse than one would hope for given that haloes are 
many-particle, extended objects which necessarily exclude 
each other on scales comparable to their radius. Correlation 
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Figure 8. The dashed curves show the bias factors between both 
the PSCz power spectrum and the dark matter power spectra in 
each simulation; the solid curves show the bias factors between 
the haloes and the dark mat ter in each simula tion. The PSCz 
bias curves are bpscz (fc) = \J -PpSCz(fe) / Pdm(^) i with a different 
Pj m for each simulation. The simulation bias curves are 6 s i m (fc) = 



haloes (&)/Fdm(fc), again running through the simulations. An 
error band arising from the errors in the power spectrum in Figure 
I13l as been put around the 64 h~ l Mpc bias curve, representative 
of the others. 



functions and power spectra appear with their theoretical 
error bars in Figs, ll2l and ll3l The correlation functions turn 
down at large scales because waves start to damp out as one 
approaches half the box size. These figures also show the 
dark matter correlation functions from each simulation; see 
Fig.|g|for a clearer view of their relationships. The error bars 
in the correlation function were calculated by splitting the 
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Figure 9. x 2 values for the four simulations as functions of 
central density cutoff p c ,min- In the 256 h~ x Mpc simulation, the 
X 2 value at p Cjm i„ = (including all detected haloes) is shifted 
to Pc,min = 0.5 so that it can appear on a log-log plot. 
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Figure 10. x 2 values for the four simulations as functions of 
halo number density. The best-fitting number densities are 0.0244, 
0.0175, 0.00873, and 0.00611 haloes per ft" 3 Mpc 3 for the 32, 64, 
128, and 256 h~ l Mpc simulations, respectively. 
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Figure 11. The sum over the smallest three simulations of x 2 as 
a function of halo number density. The best fit is at 0.0182 haloes 
per h~ 3 Mpc 3 , which corresponds to the PSCz galaxy number 
density at a depth of about 20 ft, - 1 Mpc. 



Figure 12. The best-fitting halo correlation functions (thick 
solid curves) shown individually with error bands (thin solid 
curves). The PSCz correlation function appears as a light curve 
with grey error bands, and the simulations' dark matter correla- 
tion functions appear as dotted curves. Arrows show the scales 
of the DENMAX 2 smoothing length r smoo (6/2560), and half the 
box size (6/2), where the correlation function in a box of size 6 
becomes meaningless. All softening lengths were 0.01 ft -1 Mpc. 
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Figure 13. The best-fitting halo power spectra (thick solid 
curves) shown individually with error bands (thin solid curves). 
The PSCz power spectrum appears as a light curve with grey er- 
ror bands, and the simulations' dark matter power spectra appear 
as dotted curves. 



simulation volume into octants and calculat ing correlation 
functions in each. The error at £(fc) is then ^/[Var(&(fc))]/8, 
where i runs over all octants. The same technique was used 
to calculate power spectrum error bars. 

Figure |H] shows a plot of the bias factor b(k) = 
yj -Phaioos(fc) / Pdm(k) , a measure of the difference between 
the galaxy and dark matter power spectra, for the four sim- 
ulations. The small-scale downturns are caused by the reso- 
lution of the halo-finding algorithm in each spectrum. The 
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Table 2. Fit Information. 



Box size b [h' 1 Mpc] 32 

Particle mass [lO 8 ?! -1 M Q ] 1.8 

Halo mass mj [# particles] 764 ± 322 

Halo mass ratio Wif,/wi2j> 2.8 ± 1.8 

Physical halo mass [lO 11 /! -1 M Q ] 1.4 ± 0.6, 

Circular velocity v c [km/s] 72 ± 28 

Halo number density [h 3 Mpc~ 3 ] 0.0244 

Pc,min [# particles] 122 



shape of the bias function is similar to what Kravtsov & 
Klypin (1999) found; the dependence of bias on scale must 
be similar to this to achieve a power-law form in the galaxy 
power spectrum, given the inflection in the dark matter 
power spectrum caused by the onset of non-linearity, and 
the turnover back at smaller scales due to virialization. 

Table|5|shows the central density cutoffs of the best fits, 
along with the number densities of the best-fitting popula- 
tions, and effective cutoffs in halo mass (from Fig. [1J and 
maximum circular velocity (from Fig. The quoted errors 
are the standard deviations in mass or v c of haloes with p c 
within 1/10 of a dex of p c ,min- We ignore errors arising from 
the goodness-of-fit, which we did not include since our x~ es- 
timate is not rigorous. These errors could be sizable, though, 
particularly in the 256 Mpc simulation. The minimum at 
6 particles in Fig.[§]for the 256 ft" 1 Mpc simulation's haloes 
is quite shallow. The x 2 value including all detected haloes 
(pc,min = 0) was only slightly greater than at p c ,min = 6. 
Using either cutoff, most of the haloes are right at the detec- 
tion limit (a DENMAX mass of 10 particles), so it is quite 
possible that the mass resolution in this simulation is insuf- 
ficient to pick up the truly best-fitting population of haloes. 



3.1 Discussion 

Figure |S| shows good fits to the PSCz power spectrum 
for each simulation, but the question remains: could these 
haloes from four different simulations represent the same set 
of haloes? Encouragingly, when translated into maximum 
circular velocity, the cutoff intervals mostly overlap with 
each other. The only disjoint error intervals are between the 
256 and the 32 h~ x Mpc simulations, and once again, the er- 
rors m the 256 h' 1 Mpc case are probably underestimated. 
However, we do not expect them to be exactly the same 
populations, since in a smaller box, the higher mass and 
DENMAX 2 spatial resolutions produce more haloes with 
small separations, many of which could join together in a 
lower-resolution simulation; this fact is evident in the in- 
creased small-scale range of the correlation functions in each 
simulation. Still, there are two apparent discrepancies which 
should be understood in Tabled in halo number density and 
in implied mass cutoff. 

One possible explanation for the different-looking halo 
populations is cosmic variance; i.e. that our chosen set of ini- 
tial conditions was funny in some way. For example, because 
of periodic boundary conditions, the 32 Mpc simulation 
is not big enough to contain the local supercluster from the 
PSCz cube we used to pick the initial conditions. We mea- 
sured the correlation functions of haloes in the central 16 
h~ x Mpc cube of each simulation, first applying the shifts 
and deformation tensors necessary to put the larger Simula- 
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Figure 14. Best-fitting halo correlation functions measured from 
the central cube of length 16 h~ 1 Mpc from each simulation, 
when constrained to have matching number densities and when 
individually fit. The larger halo sets have been mapped using the 
shift and deformation described at the end of i|2.1l on to the 32 
h~ 1 Mpc central region before finding the correlation function. 
To boost the number of pairs, we actually calculated the cross- 
correlation of haloes in the central boxes with all haloes in the box 
in each case. The PSCz correlation function is the dotted curve, 
and correlation functions from the 32, 64, 128, and 256 h~ 1 Mpc 
simulations appear as solid, dashed, dashdotted, and dashdotdot- 
dotted curves, respectively. We have put error bands (measured 
as previously by dividing the region into octants) around the 64 
h~ 1 Mpc correlation function, which are typical of the others. 

tions on top of the 32 h~ x Mpc simulation. Figure HH shows 
the results of this test, for both the individually-fit and num- 
ber density-constrained central density cutoffs. The central 
box is evidently undercorrelated relative to the larger boxes, 
but the correlation functions from different simulations seem 
consistent with each other, even when the sets of haloes are 
constrained to have the same number density. 

To try to understand the behavior of DENMAX mass 
and halo number density across the simulations, we iden- 
tified a few haloes by eye in the central regions of the 
four simulations. We looked for large haloes near each 
other across the simulations, and then visualized them 
with the points program written by Michael Blanton, 
at http://physics.nyu.edu/~mbl44/graphics.html Iden- 
tifying the same halo in each simulation was complicated 
by the bulk motion and distortion of the central regions, as 
discussed at the end of H2.ll We found four obvious haloes; 
the numbers of particles which comprised them in the 256 
hT x Mpc simulation were 636, 501, 100, and 89. Figure Hoi 
shows one of them ('Halo 3' from Table |3J, and how it was 
split by DENMAX 2 . While the results across simulations 
are similar, some subhaloes appear differently in different 
simulations, and some of them are questionably existent. If 
subhaloes have orbited their parent halo a few times, they 
might not be in the same place across all simulations, since 
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Table 3. Splitting of the same haloes with decreasing box size (and thus decreasing DENMAX 2 smoothing length). The halo number 
is the mass rank in the 32 h~ 1 Mpc simulation DENMAX list, and b denotes the box size, in h~ 1 Mpc. 
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Figure 15. Density profiles, as a function of distance from their 
DENMAX centres, of the four haloes identified across the simu- 
lations. The solid, dotted, dashed, and dashdot curves are from 
the 32, 64, 128, and 256 Mpc simulations, respectively. The 
halo appearing in Fig. 1161 is in the bottom-right corner. 



the simulations do quite well in tracking orbits, but not nec- 
essarily phases along them. Figure [TBI shows density profiles 
for the four haloes identified across the simulations, which 
agree quite well. In the next two sections, we will discuss 
what we learned from these haloes. 



3.1.1 Number Density 

From Table [5] the number density of best-fitting haloes 
changes by a factor approaching two between adjacent sim- 
ulations. To gauge the severity of this potential problem, 
we investigated the effect of forcing the number density of 
haloes in each box to match. Figure [TU1 shows \ 2 as a func- 
tion of halo number density. Excluding the 256 ft -1 Mpc sim- 
ulation because it was not clear that the desired set of haloes 
was well-resolved, and adding together \ 2 values from other 
simulations with equal weight, we obtained a total x 2 curve, 
which appears in Fig. 1111 The best-fitting number density 
from Fig. llll was 0.0182 h~ s Mpc 3 , which matches the num- 
ber density of PSCz galaxies that are about 20 Mpc 
from the Milky Way. Figure |7| shows the resulting power 
spectra, which are not particularly consistent, either with 
each other or with PSCz. The number density constraint 
allows smaller haloes into the 128 ft -1 Mpc set, which low- 
ers the power spectrum; it has the opposite effect on the 32 
/i -1 Mpc set, removing smaller haloes and lifting the power 
spectrum. 

However, it is not obvious that we should expect the 



same number density of haloes when resolution changes. 
Suppose we have a set of N haloes in a box of volume V. 
To measure the correlation function, we bin pairs of haloes 
by distance n to obtain P(ri), the number of pairs in bin 
i, with volume V{ri). The value of the correlation function 
in bin i is then given by: 



V P(n) 
V(n) N 2 ' 



(3) 



The RHS expresses the number of halo pairs in a bin, nor- 
malized by dividing by three quantities: the volume of the 
bin, the number of haloes, and the mean halo number den- 
sity. Now suppose that each of these haloes in fact consists 
of two subhaloes which become resolved when the resolution 
increases, and that all separations between these subhaloes 
are much smaller than the rVs. Each pair turns into four, 
and N doubles. The new correlation function £' over the 
previous range of r$'s is given by 



V 4P(rQ 
V{n) (2iV) 2 



1 + £(»■* 



(4) 



and pairs also appear in new, smaller-scale bins, imparting 
to the correlation function there a value possibly as large 
as the other £(r;)'s, depending on the bin size and distri- 
bution of subhalo distances. So in this simple case, when a 
higher spatial resolution reveals more substructure, the halo 
number density increases without changing the correlation 
function except in new, smaller-scale bins. 

Might this model resemble our simulations? Table 
[3] shows how the number of subhaloes uncovered by 
DENMAX 2 changes with box size for the four haloes we 
identified. This table includes all subhaloes, not merely the 
ones which make the central density cutoff. Also, these are 
relatively large parent haloes, more likely to have substruc- 
ture than smaller ones. This seems to be a reasonable ex- 
planation for the change in halo number density across the 
simulations. 



3.1.2 Mass cutoff 

Where nib is the DENMAX mass (in number of particles) of 
a halo in a simulation of box size b, m^lm^b should ideally 
be 8, the same factor by which the mass resolution changes. 
The values of mbjm^b for the best fits in Table [5] fall well 
short of this. For our four-halo sample, m^/m&i = 6.3±0.8, 
m64/n-i28 = 6.1 ± 0.6, and m256/wi28 = 7.3 ± 0.5, some- 
where between the ratios in Tableland the expected value 
of 8. The shrinkage of a halo with increased mass resolution 
is obvious in Fig. 1161 and comes primarily from the reduction 
in r S moo in physical units as particle mass decreases. 
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Figure 16. Comparison of a halo (the third-largest in the 32 h~ 1 Mpc simulation, with a mass of a few times 10 13 Mq) in different 
simulations, and the effect of DENMAX 2 on it. The first row shows all particles, whether they were detected by DENMAX or not, in 
a box f.8 Mpc on a side around the center of the halo, in, from left to right, the 32, 64, f28, and 256 h~ 1 Mpc simulations. The 
second row shows particles in the halo as detected by DENMAX. Note that there are many evident 'subhaloes' which do not appear in 
the second row; these were already detected as different individual haloes by DENMAX, before application of DENMAX 2 . DENMAX 2 
split the DENMAX halo into 55, 20, 3, and 1 subhalo(es) in the 32, 64, f28, and 256 h" 1 Mpc simulations, respectively. The satellite 
subhaloes appear in the third row, with circles around the subhaloes (one in each of the 128 and 32 h _1 Mpc simulations) satisfying 
the central density cutoff. The final row shows the main subhalo (also satisfying the central density cutoff) as detected by DENMAX 2 . 
Note the physical enlargement of the halo with boxsize; this is a result of the increase in the smoothing length in physical units. The 
figure was produced using Nick Gnedin's IFRIT visualization tool, at http://casa.colorado.edu/-gnedin/IFRIT/ 

4 CONCLUSIONS by imposing a central density cutoff on haloes. This cen- 

tral density cutoff implies a rough dark matter mass cutoff 
We have shown that we can reproduce the clustering prop- in PSC ^ galaxies of 10 11 " 12 M , and a cutoff in maximum 
erties of infrared-selected PSCz galaxies fairly well in simu- circular velocity of roughly 100 km/sec. Thus, it appears 
lations for scales k < 40 h Mpc -1 , near our resolution limit, that dark matter physics alone is sufficient to describe the 
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distribution of PSGz galaxies on the scales we probed. It 
is doubtful that a there is a strict central density cutoff in 
haloes which house PSCz galaxies in nature, but it seems 
that central density is a still a decent indicator of the hospi- 
tality of haloes toward nascent galaxies. While the fits look 
good, their full statistical assessment awaits a covariance 
matrix for the PSCz power spectrum. 

We have also found that the best-fitting halo popula- 
tions from the four simulations are probably consistent with 
each other. Their circular velocities are roughly consistent. 
On the other hand, their number densities and mass cutoffs 
do increase systematically with spatial resolution, but we do 
not believe that these discrepancies indicate that the popu- 
lations are necessarily inconsistent. One might in fact expect 
the halo number density to increase when higher resolution 
reveals more substructure, without affecting the correlation 
function. Also, the varying mass cutoff is at least partially 
an artifact of our halo-finding algorithm. We are currently 
considering new methods to identify haloes with fewer (per- 
haps even zero) free parameters, which we hope will coax 
the tips more confidently from their simulated dark-matter 
icebergs. 
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